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Abstract. We propose a nonconforming finite element method for iscntropic 
viscous gas flow in situations where convective effects may be neglected. We 
approximate the continuity equation by a piecewise constant discontinuous 
Galcrkin method. The velocity (momentum) equation is approximated by a 
finite clement method on div-curl form using the nonconforming Crouzeix— 
Raviart space. Our main result is that the finite element method converges to 
a weak solution. The main challenge is to demonstrate the strong convergence 
of the density approximations, which is mandatory in view of the nonlinear 
pressure function. The analysis makes use of a higher integrability estimate on 
the density approximations, an equation for the "effective viscous flux", and 
rcnormalizcd versions of the discontinuous Galcrkin method. 



Let C K , with N = 2 or 3, be a bounded polygonal domain with Lipschitz 
boundary 9f2 and let T > be a fixed final time. In this paper, we consider the 
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mixed hyperbolic-elliptic type system 

d t g + div(gu) = 0, in(0,T)xfi, (1.1) 

-fiAu- \Vdivu + Vp(g) = /, in(0,T)xf2, (1.2) 

with initial data 

£»|t=o = Qo, in ^- (1-3) 

The unknowns are the density g = g(t,x) > and the velocity u = u(t,x) € R N , 
with x £ and t G (0, T). The source term / is a given function representing 
body forces such as gravity. We denote by div and V the usual spatial divergence 
and gradient operators and by A the Laplace operator. At the boundary dfl, the 
system is supplemented with the homogenous Dirichlet condition 

u = 0, on(0,T)x9n. 

The pressure p(g) is governed by the equation of state p(g) = ag 1 , a > 0. 
Typical values of 7 ranges from a maximum of | for monoatomic gases, through 
I for diatomic gases including air, to lower values close to 1 for polyatomic gases 
at high temperatures. Throughout this paper, we will always assume that 7 > 1, 
which is the most difficult case. The viscosity coefficients /i, A are assumed to be 
constant and satisfy /1 > 0, -/VA + 2/i > 0. 

The system (1.1)— (1.2) is a gross simplification of the iscntropic compressible 
Navier-Stokes equations. It provides a reasonable approximation in situations 
where convective effects may be neglected. Solutions of (1.1)— (1.2) have also been 
utilized by Lions [12] to construct solutions of the isentropic compressible Navier- 
Stokes equations. Regarding the mathematical theory, the semi-stationary system 
(1. 1)— (1.3) has been analyzed by Lions [12, Section 8.2], among many others. More 
precisely, he proves the existence of weak solutions and provide some uniqueness 
and higher regularity results. 

In the literature one can find a variety of numerical methods for the compressible 
Stokes and Navier-Stokes equations. However, there are few results with reference 
to the convergence properties of these methods, especially in several dimensions. 
In one dimension, we refer to the works of Hoff and his collaborators [15, 16, 
17]. These results apply to the compressible Navier-Stokes equations written in 
Lagrangian form and requires the initial density to be of bounded variation. In 
several dimensions there arc a few very recent results. In [7, 8], the authors present 
a convergent finite element method for a Stokes model. This model is a stationary 
version of (1.1)— (1.2). In their finite element method the approximation spaces for 
the density and velocity are the same. Moreover, their method is based on the 
standard weak formulation of the velocity equation (1.2). Since the finite element 
space is non-conforming, this approach may not preserve the div-curl structure of 
the continuous system. This complicates the convergence proof. In [7, 8], additional 
stabilization terms are needed in the discretization of the continuity equation (1.1). 
In [11], we construct a convergent mixed finite clement method for (1.1)— (1.2). 
However, this method is based on a vorticity formulation of the velocity equation, 
which is only valid for the Navicr slip boundary condition: 

u ■ v = 0, curl it x v = 0, on d£l. 

In addition, the velocity is approximated by a ff(div) (Nedelec) element. 

We now outline the numerical method proposed in this paper. First of all, 
the density g is approximated by piecewise constants in the spatial and temporal 
variables. For the approximation of the velocity u we utilize the Crouzeix-Raviart 
element space [4] in the spatial variable, denoted by Vh (fi), and piecewise constants 
in the temporal variable. Hence, the numerical method is nonconforming in the 
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sense that Vh <f. Wq' 2 (Q,). In what follows, we mostly suppress the time variable t 
and refer to subsequent sections for precise statements. For the continuity equation 
(1.1) we make use of a discontinuous Galerkin method. To achieve stability, the 
numerical fluxes are evaluated in the upwind direction dictated by the velocity. 
However, since the velocity space is not continuous across element faces, average 
velocities are used in this discretization. Our discontinuous Galerkin method is 
equivalent to a standard finite volume method for the continuity equation [6, 9]. 
In [11], we use a similar discontinuous Galerkin method with the velocity in the 
div conforming Ncdclec space of the first order and kind. Since the method used 
herein only depends on the average normal velocity at faces, the approximations 
constructed by this method are also solutions to the discrete continuity equation 
of [11]. More precisely, if the pair (gh,Uh) solves the discrete continuity equation 
proposed herein, then (gh,^h u h) is a solution to the discrete continuity equation 
of [11], where 11^ is the canonical interpolation operator onto the div conforming 
Nedelec space of first order and kind. As a consequence, several of the favorable 
properties of the method in [11] continue to hold for the continuity method herein. 
In particular, renormalized formulations, weak time-continuity, and consistency 
bounds are readily obtained by exploiting this connection. 

To discretize the velocity equation (1.2) we bring into service a non-standard 
finite clement formulation, which starts off from the identity 

/ DuDv dx — / curl u curl v -\- div u div v dx, (1.4) 
Jn Jn 

valid for all u G W 1,2 (Q). However, since the velocity space is nonconforming, this 
identity does not hold discretely, but we insist on utilizing the right-hand side of 
(1.4) as a starting point for discretizing the velocity equation. Utilizing the form 
on the right-hand side, it is possible to split the curl part of the Laplacian away 
from the divergence part. By setting v = Vs, we obtain the divergence part, while 
v = curl 77 gives the curl part. Of course, to satisfy boundary conditions, this 
argument must be localized. Discretely, this still holds for the element space Vh 
since this admits the exact orthogonal Hodge decomposition 

V h = curia + VS h . 

Hence, the curl and divergence part of the Laplace operator can be separated by 
using test functions Vh = curlOi, Ch <= Wh or Vh = Vs/j, Sh € Sh- This property 
lies at the heart of the matter in the upcoming convergence analysis. 

Contrasting with the standard situation in which the left-hand side of (1.4) is 
used, a discretization based on the right-hand side of (1.4) does not converge unless 
additional terms controlling the discontinuities of the velocity are added, cf. Brenner 
[2]. The standard discretization of the Laplacian (based on the left-hand side of 
(1.4)) leads to a L 2 bound on VhUh, where is the gradient restricted to each 
element E. For the velocity space Vh, this bound actually controls the jump of Uh 
across faces. This in turn, is sufficient to conclude that V^m/i — 51 Vu as h — > 0. 
When discretizing the Laplacian based on the right-hand side of (1.4), one obtains 
L 2 bounds on curl/j Uh and div/jit^, where curU and div/j denotes the curl and 
divergence operators, respectively, restricted to each clement E. The jump of Uh 
across faces is not controlled by these terms. In fact, Vh contains non-zero functions 
for which both div^ and curl^ are zero. For this reason, extra terms controlling the 
jump of Uh across faces need to be added. 

In choosing these terms we are inspired by the work of Brenner [2] , which deals 
with two-dimensional elliptic operators of the form "curl curl — /3V div" . To be more 
precise, our finite element method for the velocity equation (1.2) seeks Uh € Vh(fl) 
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f h v h dx, Vv h e V h (Cl), 



such that 

/ /ucurl/j «/, curl/, u h + [(^ + A) div/, u h - p(gh)} div h v h dx 
Jn 

+ M / I"' 1 ' H ' Hlr + x ^lr H x "lr dS ( x ) (1.5) 

-/ 

Jn 

for some fixed e S (0,1), where ph and fh are given piecewise functions on with 
respect to a tetrahedral mesh Eh with elements £7. Moreover, denote the set of 
internal faces, and |-] r denotes the jump across a face V E T^. The scaling factor 
h e is required to prove convergence of the finite element method. Of course, the 
size of e will affect the accuracy of the method [2] and should be fixed very small 
in practical computations. 

For any fixed h > 0, let (Qh,Uh) = (gh,Uh)(t, x) denote the numerical solution to 
the compressible Stokes system. Our goal is to prove that {(Qh,Uh)} h >o converges 
along a subsequence to a weak solution. The main challenge is to show that the 
density approximations Qh, which a priori is only weakly compact in I? , in fact 
converges strongly. Strong convergence is needed when sending h — > in the 
nonlinear pressure function. It is this issue that motivates the above nonconforming 
finite element method. Since the finite element space Vh is piecewise linear and 
totally determined by its value at the faces, Green's theorem yield 

div,, U\v = div v, curl/, Ii\v = 11^ curl v, 

where 11^ is the canonical interpolation operator onto Vh and 11^ is the L 2 projec- 
tion onto piecewise constants. Consequently, the projection of a divergence or curl 
free function is again (piecewise) divergence or curl free. Using this, we see that 
the function Vh — II^VA _1 p^ is a solution to the div-curl problem 

div h v h = Qh, curl/, v h = 0, 

away from the boundary. By using Vh as test function in (1.5), the curl term 
vanishes, while the remaining terms constitute the so-called effective viscous flux 
P e s{Qhi u h) = p(ffh) ~ (A + fi) div Uh, the source term, and the jump terms. The 
latter terms are shown to converge to zero. Using this, we are able to prove following 
weak continuity property: 

lirn J J P e s{eh,Uh) Qh4> dxdt = J J P e g g<j) dxdt (P e s, Q are weak L 2 limits), 

(1-6) 

for all 4> € Co° ■ This is the main ingredient in the strong convergence proof for the 
density approximations Qh- The argument is inspired by the work of Lions on the 
compressible Navier-Stokes equations, cf. [12]. 

If we instead of (1.5), discretize the Laplacian based on the left-hand side of 
(1.4), then the above analysis becomes more involved. In particular, it seems diffi- 
cult to establish the key property (1.6). In this case, we would need to establish 



J jVhUhVhK [V^ Qh ]-^ h u heh dxdt 



0, as h — > 0, 



which is intricate since all the involved quantities are only weakly convergent. 

The remaining part of this paper is organized as follows: In Section 2, we first 
introduce some relevant notation and state a few basic results from analysis. Next, 
we formulate our notion of a weak solution. Finally, we introduce the finite element 
spaces and derive some of their basic properties. In Section 3, we present the 
numerical method and state our main convergence result. This section also provides 
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a result regarding the existence of solutions to the discrete equations. In Section 4, 
we derive stability and higher intcgrability results. Section 5 is devoted to proving 
the convergence result stated in Section 3. 

2. Preliminary material 

2.1. Functional spaces and analysis results. We denote the spatial divergence 
and curl operators by div and curl, respectively. As usual in the two dimensions, we 
denote both the rotation operator taking scalars into vectors and the curl operator 
taking vectors into scalars by curl. 
We will make use of the spaces 

w dw < 2 (n) = {v e L 2 (n) : diw e L 2 (n)} , 

W cu * h2 (n) = {ve L 2 (n) : cuvlv e L 2 (fl)} , 

where v denotes the unit outward pointing normal vector on dfl. If v € W dlv -' 2 (Q) 
satisfies v ■ v\qq = 0, we write v e Wq 1V ' 2 (Q). Similarly, v e WQ Url ' 2 (fi) means 
v € W curl ' 2 (fi) and v x v\ dn = 0. From [10], 

w 1,2 (fi) = w cur1 - 2 n w div ' 2 . 

The next lemma lists some basic results from functional analysis to be used in 
subsequent arguments (for proofs, see, e.g., [5]). Throughout the paper we use 
overbars to denote weak limits, in spaces that should be clear from the context. 

Lemma 2.1. Let O be a bounded and open subset of M M with M > 1. Suppose 
g: R — > (—00,00] is a lower semicontinuous convex function and is a 

sequence of functions on O for which v n — ^ v in i 1 (0), g(v n ) € ^(O) for each 
n, g(v n ) — 1 g(v) in i 1 (0). Then g(v) < g(v) a.e. on O, g(v) £ ^(O), and 
J g{v) dy < liminfn^oo J g(y n ) dy. If in addition, g is strictly convex on an 
open interval (a,b) C K and g{v) — g(v) a.e. on O, then, passing to a subsequence 
if necessary, v n (y) -> v(y) for a.e. y G {y e O \ v(y) e (a, b)}. 

Let X be a Banach space and X* its dual. The space X* equipped with the 
weak-* topology is denoted by X* eak , while X equipped with the weak topology 
is denoted by X weak . By the Banach- Alaoglu theorem, bounded balls in X* are 
o~(X*, X)-compact. If X separable, the weak-* topology is metrizable on bounded 
sets in X* , which makes it possible to consider the metric space C ([0, T]; X* cak ) 
of functions v : [0, T] — > X* that are continuous with respect to the weak topology. 
We have v n -» v in C ([0, T]; X* eak ) if (v n (t),<f>) x *,x -» (v(t),<f>)x* ,x uniformly 
with respect to t, for any (f> e X. The succeding lemma is a consequence of the 
Arzela-Ascoli theorem: 

Lemma 2.2. Let X be a separable Banach space, and suppose v n : [0,T] — > X* , 
n = 1,2,..., is a sequence for which ||w„|| loo q T y x *) — ^, for some constant C 
independent of n. Suppose the sequence [0, T] 3fn (v n (t), &)x*,x, n — 1,2,..., 
is equi- continuous for every $ that belongs to a dense subset of X . Then v n belongs 
to C ([0,T];X* eak ) for every n, and there exists a function v £ C ([0, T]; A* eak ) 
such that along a subsequence as n — > 00 there holds v n — > u in C ([0, T]; A* oak ). 

Later we frequently obtain a priori estimates for a sequence {v n } n>1 that we 
make known as "v n Gb A" for a given functional space X. What this really means 
is that we have a bound on ||i> n ||x that is independent of n. 

The following discrete version of a lemma due to Lions [12, Lemma 5.1] will 
prove useful in the convergence analysis. 

Lemma 2.3. Given T > and a small number h > 0, write (0, T] = U^l 1 (ffe_i, tk] 
with tk = hk and Mh = T. Let {fh}'h > Q, {ffh}h>o ^ e ^ wo sequences such that: 



6 



K. H. KARLSEN AND T. K. KARPER 



(1) the mappings t i— > gh(t,x) and t i— > fh(t,x) are constant on each interval 
(tfc-i,tfc], k = 1,...,M. 

(2) {fh}h>o and {9h}h>o converge weakly to f and g in L Pl (0, T; L 91 (f2)) and 
L P2 (0, T; L 92 (f2)), respectively, where 1 <Pi,ffi < oo and 



1 1 

— + — 

Pi 



1 1 

— + — 
91 92 



1. 



(3) £/ie discrete time derivative satisfies 

e b L\0,T;W 



g h {t,x)-g h {t-h,x) ^ rl , n ^ n7 . u 



(fi)) 



(4) ||//i(*,a;) - fh(t,x- OIU"2(o,T;L92(n)) -» as |£| -> 0, uniformly in h. 
Then ghfh — ^ 5./ in sense of distributions on (0,T) x fL 

Proof. Let us introduce an auxiliary piecewise linear function by setting 

5fc (V) = 5/i(*fe) + h^it-tk) (gh(tk+i) - 9h{tk)) , t e (ife,ife+i], 
for fc = 0, . . . , M — 1. Using property (3), 



gh{t,x) -g h (t-h,x) 
h 



dxdt 



(2.1) 



< Ch\\4>\\ L oo {0!T . W i,oc (n)) , (f> e C °°(fi). 

Thus, (g/j — .Oh) — ^ as in the sense of distributions on (0, T) x £1 as h — > 0. 
Next, we write 

3/i //» = 5/i//i + (5/i ~ <7h)/h- 
By requirement (3), d t gh ^(OjT; W _1,1 (fi)). This and requirement (4) allow 
us to apply a lemma due to Lions [12, Lemma 5.1], yielding 

fh9h fg, 

in the sense of distributions on (0, T) x ft as h — > 0. 

It only remains to prove that (gh — gh)fh ^ in the sense of distributions. 
For this purpose, set = fh * K e , where k £ is a standard smoothing kernel and * 
denotes the convolution product. We write 

(gh - 9h)h = {gh - 9h)fh + {gh - 9h)(h - ft)- 

Now, requirement (4) yields 

Wfh - fh\\LP2(0,T;Li2(Q)) — > as 6 — > 0, 

uniformly in /i, and hence 

lim lim / (gh - 9h)(fh - fh)<t> dxdt = °- 
Thus, the proof is complete provided that 

lim lim / (g h - 9h)fh^> dxdt = °- 
e^Oh^oj J n 

By a calculation similar to (2.1) we see that 

cT 



(9h - 9h)fh dxdt 



< ft "2 C||//[||£P2(0,T;W' 1 .°°(fi))) 



where we have also applied Lemma 2.10 (below) to the time variable. From this 
we can conclude that (gh — 9h)fh ^ in the sense of distributions as h — > 0. This 
brings the proof to an end. □ 
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2.2. Weak and renormalized solutions. 

Definition 2.4 (Weak solutions). A pair of functions (g, u) constitutes a weak 
solution of the semi-stationary compressible Stokes system (1.1)— (1.2) with initial 
data (1.3) provided that: 

(1) {g,u) e Z/°°(0, T; £ 7 (f2)) x L 2 (0, T; W 1,2 (fl)), 

(2) d t g + div(gu) = in the weak sense, i.e, V</> € C°°([0,T) x U), 

[ [ g{4> t + uD<j)) dxdt+ [ g <t>\ t=a dx = 0; (2.2) 
Jo Jn Jn 

(3) —fj,Au — \D div u + Dp(g) = f in the weak sense, i.e, V0 e C o °°([0,T)xfi), 

[ [ p\7u\7(f> + [(n + Xdiv u-p(g)} div <f> dxdt = [ f f<j>dxdt. (2.3) 
Jo Jn Jo Jn 

For the convergence analysis we shall also need the DiPerna-Lions concept of 
renormalized solutions of the continuity equation. 

Definition 2.5 (Renormalized solutions). Given u E i 2 (0, T; W 1,2 (fi)), we say 
that g € L°°(0,T; L 7 (n)) is a renormalized solution of (1.1) provided 

B(g) t + div (B(g)u) = b(g) div u in the sense of distributions on [0, T) x Q, 

for any B G C[0, oo) fl C 1 (0, oo) with B(0) = and b(g) := B'{g)g-B{g). 

We shall need the following well-known lemma [12] stating that square-integrable 
weak solutions g arc also renormalized solutions. 

Lemma 2.6. Suppose (g,u) is a weak solution according to Definition 2.4-. If 
g £ L 2 ((0,T) x ft)), then g is a renormalized solution according to Definition 2.5. 

Remark 2.7. Regarding the continuity equation and the definitions of weak and 
renormalized solutions, we are requiring the equation to hold up to the boundary. 

2.3. On the equation div v = /. Solutions to the following problem are vital to 
the upcoming convergence analysis: 

divu = / in fl, v = on dfl. (2.4) 

If / e L P (Q) with J Q f dx = 0, then a solution to (2.4) can be constructed through 
the Hodge decomposition 

v = Vs + curl£, 

where s e H 2 (Q) solves the Neumann Laplace problem, i.e., 

As = / in CI, Vs ■ v = on dCt, 

and £ G H 2 (Q) is determined such that v\qq = (cf. [1]). Such a solution can be 
constructed using the Bogovskii solution operator [5]. Here, we define the solution 
operator B [■] : Lg(fi) -» Wo' p (n) as one of the solutions to the problem 

dWB [(j>\ = 4> in fl, B[4>] = on dfl. (2.5) 

We shall need solutions v satisfying cwclv = 0. Clearly, this is not compatible 
with the Dirichlet boundary condition. However, locally curl free solutions can be 
constructed using the operator A[-] : L p (Cl) — > W 1 ^^), 

A[<j>] = VA- 1 [<j>], (2.6) 

where A -1 is the inverse Neumann Laplace operator. 
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2.4. Finite element spaces and some basic properties. Let Eh denote a shape 
regular tetrahedral mesh of ft. Let T h = {r e Th ■ T <£_ dfl} denote the set of 
internal faces in Eh- We will approximate the density in the space of piecewise 
constants on E h and denote this space by Q^(fi). For the approximation of the 
velocity we use the Crouzeix-Raviart element space [4]: 

Vfc(fi) - |« fc g L 2 (n) : v h \ E e P?{E), VE g E h , H] r dS(x) = o, vr e r£ j , 

(2.7) 

where J-] r denotes the jump across a face T. To incorporate the boundary condition, 
we let the degrees of freedom of Vh(fl) vanish at the boundary. Consequently, the 
finite clement method is nonconforming in the sense that the velocity approximation 
space is not a subspacc of the corresponding continuous space, Wq' 2 (H). 
We introduce the canonical interpolation operators 

Hi : W^ 2 (n) - Vfc(fi), Il£ : L 2 (n) -> Q h (fi), 



defined by 



J uXv h dS(x) = J Vh dS(x), vr G r h , 
/ dx= I <j) dx, \/e e £ h . 



(2i 



Then, by virtue of (2.8) and Stokes' theorem, 

div^n^t; = n^divu, curl*, = 11% cmlv, (2.9) 



for all v e W 1,2 (il). Here, curl^ and divh denote the curl and divergence operators, 
respectively, taken inside each element. 
Now, (2.9) immediately gives 

div h UlB [q h ] = q h , Mq h G Q/,(fi) H j jT % rfa; = J , 

and, away from the boundary, 

div h III A [q h ] = q h , curl/, Iljf.4 = 0, Vg fc G Qh(ft), 

where B [■] and *4 [•] are defined in (2.5) and (2.6), respectively. Consequently, this 
configuration of elements enables us to construct discrete analogs of the continuous 
operators (2.5) and (2.6). 

We associate to the space V/j(f2) the following semi norm and norm: 

\vh\v h = \\carl h v h \\ 2 L 2 (n) + \\ div h v h \\ 2 L2{n) 

+ h " X (H l Vh ■ v 1t \\1-{T) + II H x i/] r ||| 2(r) ) , (210) 
rer h 

\\vh\\v h = \\Vh\\h(0) + \ v h\v h - 

Let us now state some basic properties of the finite element spaces. We start by 
recalling from [3, 4] a few interpolation error estimates. 

Lemma 2.8. There exists a constant C > 0, depending only on the shape regularity 
of E h and such that for any 1 < p < oo, 

||n^-0|U P(n) <Ch\\V<P\\ LP(Q) , 

||nj> - w||iP(n) + h\\V h (ILXv - v)\\ LP( a) < ch s \\V s v\\ LP{E) , s = 1,2, 

for all <j> e W 1,p (£i) and v G W S ' P (E). Here, \7 h is the gradient operator taken 
inside each element. 
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By scaling arguments, the trace theorem, and the Poincare inequality, we obtain 

Lemma 2.9. For any E G Eh and <f> G W ' (E), we have the following inequalities: 

(1) trace inequality, 

UWmv) < ch E h (UWlhe) + h E \\V<]>\\v{E)) , vr e r h n dE. 

(2) Poincare inequality, 

'~w\L 



dx 



L 2 (E) 



<Ch E \\V<P\\ L 2 {E) . 



In both estimates, h E is the diameter of the element E. 

Lemma 2.10. There exists a positive constant C, depending only on the shape 
regularity of Eh, such that for 1 < q,p < oo and r = 0, 1, 

Uh\\ W r, P{E) < Ch- r + ^f-f > Uh\\ Lq{E) , 

for any E G Eh and all polynomial functions 4>h G Fk(E), k = 0, 1, . . .. 

Lemma 2.11. Let {vh}h>a be a sequence in Vh(£l). Assume that there is a constant 
C > 0, independent of h, such that \\vh\\v h < C. Then there exists a function 
v G WQ 1,2 (f2) such that, by passing to a subsequence as h — > if necessary, 

Vh — 1 v in L 2 (Q), cwclhVh — 1 curl I? in L 2 (Q), div/j Vh drvv in L 2 (Q). 

Proof. As ||w/i||v h i s bounded independently of h, it follows that Vh Gh L 2 (fl), 
curU Vh &b L 2 (ty, and div/j Vh G& L 2 (f2). Thus, we have the existence of functions 
v G L 2 (fl), £ G i 2 (^), and £ G L 2 (f2) such that, by passing to a subsequence if 
necessary, 

v h v, curlft v h £, div ft C 

Once we make the identifications £ = curlv and £ = divu, the proof is complete. 
Fix any G Wj' 2 (fi). An application of Green's theorem yields 



curl ft Vh<t> dx = 

EGE h 



Vh curl </> da; 



as 



x z/) dS(x) 



= / Vh curl da; + V, / ^[v/,X!/] r dS(x). 
Jn rerI Jr 



By sending h — ► in the above identity, we discover 



/ ^4>dx= i v curl rfa; + lim / ^[^xj/L dS(a;). (2.11) 



Utilizing the bound 



rer' ^ 



(2.12) 



cf. (2.10), and the second condition in (2.7), we control the last term of (2.11): 



V f 4>lv h xv} T dS( 



x) 



rerf 



6 r ) \v h x v\ T dS(x) 
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< ft - 5 



[Vh ~xv\ T 



£ ft j^-M 2 dS(x) 



where {</>r}rer h is a given set of real numbers. For each T e T^, let us take 
4>y '■= J E <fi dx, where E is arbitrarily fixed as one of the two elements sharing 
the edge T. Now, using Lemma 2.9 and (2.12), we deduce 



V /^xi/] r dS(x) 
rer< J * 



< Ch-ih\\W<t)\\ L 2 



(SI)- 



Hence, 



lim 



V ( ct>iv h xv\ T dS(x) 



= 0. 



By (2.11), this shows that 



/ dx = / v 

J si J si 



curl 6 dx, 



and so our claim follows, i.e., £ = curl v. 

By almost identical arguments we find that £ = divu. 



□ 



The following lemma provides us with an estimate of the blow-up rate of VhVh, 
for any element Vh S Vh(fi). 

Lemma 2.12. There exists a positive constant C, depending only on the shape 
regularity of Eh and the size of ft, such that 



rllz, 2 (r) 



for all v h e V h (Q). 



Proof. By the linearity of v h \ E , Av h \ E = \/E e E h . Using this we can apply 
Green's theorem to deduce the bound 



\\^hV h \\ 2 L 2 {n) = £ / ^hVh-^hVh dx= ^2 / VJv h ■ v)v h dS(x) 
= £ / (VK] r -f)«fc ^(z) 

< £ / IVKlr-HKI dS(x)=:I. 

To obtain the third equality we have used that the average of is continuous 
across internal faces. Since Vh € Vh(Cl), we know that V[w/,] r is constant for all 
internal faces T £ T^. Moreover, there must exist a point br € I\ for every rer[, 
such that [vh(&r)]r = 0- By this an d the Cauchy-Schwartz inequality, we deduce 

I<C £ ^KII^(r)IIM r ll^(r) 



^Cft-^ill^Hi^n) £ ft 6 - 1 !!^] 



r lli 2 (r) 
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The last inequality is achieved thanks to the trace inequality (1) in Lemma 2.9, 
together with Lemma 2.10. □ 

Using the previous lemma, we can now establish a Poincare inequality and a 
spatial compactness estimate. 

Lemma 2.13. There exists a positive constant C , depending only on the shape 
regularity of Eh and the size of Q,, such that for any ^eM 2 

\\v h (-)-v h (--0\\LHn ( )<C\&-*\v h \ Vh{Q) , Vv h eV h (fl), (2.13) 

where = {x e ft : dist(a;, dCl) > £}. Moreover, 

IKMfi) < CKk(fi), Vv h EV h (n). (2.14) 

Proof. Fix h > 0, and select an arbitrary function Vh in Vh(fl). For any £ we 
construct a new mesh Gh such that each G e Gh is a subset of one and only one 
element E e (e.g., we can divide each element of E € -E^ into a number of 
smaller elements). Moreover, we construct this new mesh Gh such that 

C" 1 ^ < h G < C\£\, VGeGh, (2.15) 

where he denotes the diameter of the new element G and the constant C depends 
only on the shape-regularity of Eh- 

Now, let V[ £| (f2) denote the Crouzeix-Raviart element space on Gh and denote 
by n^| : Vh(fl) — > Vj ? |(n) the canonical interpolation operator associated with 

Denote by fti^i the maximal element diameter in G/j. From standard properties 
of the Crouzeix-Raviart element [14], we have 

\\U^v h (x)-U^v h (x- 0lll 2( n,) <(4l + KI 2 ) E H Vn m^H^(G) 

<(4i + iei 2 )l|v^||i 2(n) , 

where the second inequality follows from the properties of the operator II 7£, . 

Lemma 2.12 and the bounds (2.15) allow us to conclude the following estimate: 

\\U^Vh(x) - U^v h (x - Olli^nt) 

<h^-Hh^ + \e)\\vh\\ LHn) I E ^-'HMrlll^r) 

<^| 1_i ll^ll^(n) ( E fe£_1 IIKlrlli»(r) 

Keeping in mind that vh\e € W 1,2 (.E), V.E G i?/,, we can apply Lemma 2.8 and 
the previous estimate to obtain 

K(-)-«fc(--0lli* ( n t ) 

< 2|K - n^|W fc ||i 2{ „ e) + ||n^0r) - n^^^ - Olli^) 

(2.16) 

<C\H\ l -Hv h \\ L , (n) \ ]T ^IIHlrllW/ 
vrer£ 

Next, denote by v h * % the extension of Vh by zero to all of R N . By the previous 
calculations, we conclude that vff 1 satisfies (2.16) with fi^ replaced by 1* (keep 
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in mind that the jump terms are only summed over internal faces). Thus, we can 
fix |£| large in (2.16) to discover that 



I cxt||2 
\\ v h \\L 2 (n) 



and hence 



Klli'(n) < Cldiam^)! 1 "* ]T ^"'H Wrll^(r) I < C\ diarn^)! 1 "* \v h \ Vh 

Ver' h J 

which is (2.14). 

Finally, setting (2.14) into (2.16) gives (2.13). □ 
We end this section with 

Lemma 2.14. There exists a constant C > 0, which depends only on the shape 
regularity of Eh and the size of ft, such that for any Vh € Vh(fl) 

< C/i5||t7 h || v . h(fl) ||Vw||i2 (n ), Vto e Wo' 2 (f2). 
Proof. Using the Holder inequality, 

]T V- 1 [ {v h -y\ T lnXwv] r dS(x) 

< /it I V tf" 1 f \v h ■ vf T dS{x) 

x I h ^ I \ U hW-w\ 2 dS(x)\ . 
\EeE h / 

To obtain the last inequality, we have applied the calculation 

KHr = K n J»U+ -«> + «,- (nj»| ejs - | 2 

< |(il»| 9e+ - ™| 2 + |(nJ»| aB - - w\\ 

where E + and E~ are the two element sharing the face T. 
By using (1) in Lemma 2.9, we further deduce that 



rer h 



< 



Vevl J * 



dS(ar) 
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x £ h- 2 \\Ulw-w\\l 2{E) + \MUlw-w)\\l 2{E) 

\E£E h / 

< hi C\\v h \\ Vh \\Vw\\ L 2( n) , 

where the last inequality follows from Lemma 2.8. 
By analogous calculations for the tangential jumps, 



£ K-- 1 [ lv h x uj T \n v h w x u\ T dS(x) 
rerl Jr 



< C\\v h \\ Vh \\Vw\\ L 2 m . 



This concludes the proof. 



□ 



3. Numerical method and main result 

Given a time step At > 0, we discretize the time interval [0,T] in terms of the 
points t m — mAt, m = 0, . . . , M, where it is assumed that MAt — T. Regarding 
the spatial discretization, we let {Eh}h be a shape regular family of tetrahedral 
meshes of il, where h is the maximal diameter. It will be a standing assumption 
that h and At are related like At = ch for some constant c. By shape regular we 
mean the existence of a constant k > such that every E e E h contains a ball of 
radius > — , where \ie is the diameter of E. Furthermore, we let Th denote 
the set of faces in Eh. Throughout the paper, we will use "three dimensional" 
terminology (tetrahedron, face, etc.) when referring to both the three dimensional 
case and the two dimensional case (triangle, edge, etc). 

On each element E e Eh, we denote by Q(E) the constants on E. The functions 
that are piecewise constant with respect to the elements of a mesh Eh are denoted 
by Qh(Q)- We denote by Vh(Q) the Crouzeix-Raviart finite element space (2.7) 
formed on Eh- To incorporate the boundary condition, we let the degrees of freedom 
of Vh(Q vanish at the boundary: 

Jv h dS(x) = o, yrer h ndn, Vv h ev h {fi). 

We shall need to introduce some additional notation related to the discontinuous 
Galerkin method. Concerning the boundary dE of an clement E, we write / + for 
the trace of the function / achieved from within the element E and /_ for the trace 
of / achieved from outside E. Concerning a face T that is shared between two 
elements E_ and E + , we will write /+ for the trace of / achieved from within E + 
and /_ for the trace of / achieved from within E_ . Here and E + are defined 
such that v points from to E + , where v is fixed (throughout) as one of the 
two possible normal components on each face T. We also write |/] r = /+ — /_ for 
the jump of / across the face T, while forward time-differencing of / is denoted by 
|ymj = pn+i _ j m and gh jm = IO. Tne S ct of inner faces of T h will be denoted 

b y r£ = {rei\ ; r£c>n}. 

Definition 3.1 (Numerical scheme). Let {&h( x )} h>o ^ e a sequence (of piecewise 
constant functions) in Qh(ty that satisfies p° > for each fixed h > and g° h — > go 
a.c. in n and in L 1 ^) as h -> 0. Set f h (t, •) = /£"(•) := ^ J^-i n^/(s, •) ds, for 
t e (i m _i,i m ), m=l,...,M. 
Determine functions 



eQ fc (fi)xV h (fi), m=l,...,M, 
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such that for all <f>h G Qh(Q), 

d?( g Z)d> h dx = AtJ2 I (? m K'^ + W'4)Wr*. (3.1) 



and for all Vh & Vh(fl), 

/ n curl/, u™ curl/j v h + [(fi + A) div h u™ - p(Qh)] div/, v h dx 
Jn 



I 



+ » E ^ / K* • "lr K ' Hlr + K* x "lr K * i/] r dS(x) (3 2) 
K l v h dx, 

in 

for m = 1, . . • , M . 

In (3.1), we have introduced the notation 

{u h -u)i=QrJu h -udS{x)j , (3.3) 
where a + — max(a, 0) and a~ = min(a, 0). 

Remark 3.2. Since the normal velocity components (w™-^) are discontinuous across 
element faces, the continuity method (3.1) approximates (gu ■ v) using instead the 
average normal velocity ^ J r ( M ™ ' v ) dS(x), cf. (3.3), and traces of g are taken in 
the upwind direction with respect to the average normal velocity. 

We now make an observation that will simplify the subsequent analysis. Let 
Afh(ty denote the lowest order div conforming Nedelec finite element space of the 
first kind [13, 11] on E h . In two dimensions, A/"h(fi) is the Raviart-Thomas space. 

We will need the interpolation operator 1T^ : Vh(Q) — > defined by 

J (i# v h ) ■ v ds(x) = J v h -v dS{x), vr g IV 

Then, by definition, the interpolated velocity 

(3-4) 

satisfies 

= (<- v ) h = «■")*. (3-5) 

where the first equality is valid since • v is constant on each r e IV Since both 
element spaces have piecewise constant divergence, a direct calculation yields 

divw™ = di\ h v%. (3.6) 

Now, setting (3.5) into the continuity method (3.1) leads to the relation 

d?( Q %)cj> h dx = At I ^(^•") + + ?K"')")Wr dS(x), (3?) 

for all <f>h € Hence, we can think of the pair (g™,it™) as a solution to a 

continuity method in which is used to approximate the velocity. In fact, 

(3.7) is the method examined in [11]. We will frequently utilize (3.7), instead of 
(3.1), to easily obtain properties of our continuity approximations. 

For each fixed h > 0, the numerical solution {{q™ ' u h')}m=o * s extended to the 
whole of (0, T] x il by setting 

{g h ,u h ){t) = {g k a ,u k n ), te(V!,g, m=l,...,M. (3.8) 

In addition, we set £>/j(0) = 
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3.1. Main result. Our main result is that, passing if necessary to a subsequence, 
{(Qh, u h)}h>o converges to a weak solution. More precisely, we will prove 

Theorem 3.3 (Convergence). Suppose f e L 2 ((0,T) x SI), and g E L^(Sl) with 
7 > 1. Let {{gh, u h)}h>o be a sequence of numerical solutions constructed according 
to (3.8) and Definition 3.1. Then, passing if necessary to a subsequence as h — ► 0, 
u h — »■ u in L 2 (0,T; L 2 (Sl)), QhUh QU in the sense of distributions on (0,T) x Si, 
and Qh — > Q a.e. in (0,T) x SI, where the limit pair (g,u) is a weak solution as 
stated in Definition 2-4- 

This theorem will be a consequence of the results proved in Sections 4 and 5. 

3.2. The numerical method is well— defined. We now turn to the existence of 
a solution to the discrete problem. However, we commence with the following easy 
lemma providing a positive lower bound for the density. 

Lemma 3.4. Fix any m = 1,...,M and suppose g™^ 1 € Qh{Sl), u™ e Vh(Sl) 
are given bounded functions. Then the solution g™ <G Qh(Sl) of the discontinuous 
Galerkin method (3.1) satisfies 

minol n > minpT 1-1 [ t . . ~ n ) . 

xm^ h ~ xen eh \l + At\\div h uf\\ L oo {n) J 

Consequently, if g"^ 1 ^) > 0, then g^(-) > 0. 

Proof. Let u™ be given by (3.4). Then, since (p™,«™) satisfies (3.7), Lemma 4.1 
in [11] can be applied. This concludes the proof. □ 

Lemma 3.5. For each fixed h > 0, there exists a solution 

( Q Z,u™)eQ h (Sl)xV h (Sl), e n-)>0, m=l,...,M, 
to the discrete problem posed in Definition 3.1. 

Proof. As in the proof of [11, Lemma 4.2], the existence of a solution is established 
using a topological degree argument. By the arguments corresponding to those 
in [11, Lemma 4.2], one reduces the problem to proving existence of a solution 
Uh € Vh(Sl) to the linear system: 

a(u h , v h ) := / ^curlft u h curU v h + (m + A) div^ u h div ?l v h dx 
Jn 

+ n ^2 fee_1 / l u h ■ v lv l v h ■ ^lr + l u h x ^lr l v h x "Ir dS ( x ) 

= / gv h dx, \/v h e V h (Sl), 
Jn 

(3.9) 

where g £ L 2 {Sl) is given. 

Now, the bilinear form a(-,-) is clearly bounded on the space Vh(Sl) x Vh(Sl) 
equipped with the norm || • ||v h - By an application of the Poincare inequality 
(2.14), we also have the existence of a constant C, independent of h, such that 

a(u h ,u h ) > C\\u h \\lr h . 

Hence, the bilinear form a(-,-) is coercive on Vh, and the existence of a function 
Uh € Vh(Si) satisfying (3.9) follows. 

Since the remaining part of the topological degree argument is very similar to 
that found in [11, Lemma 4.2], we omit the details. □ 
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4. Basic estimates 

In this section we establish various a priori estimates for the discrete problem 
given in Definition 3.1, including a basic energy estimate and a higher integrability 
estimate for the density approximations. 

We begin with a renormalizcd formulation of the continuity method (3.1). 

Lemma 4.1 (Rcnormalized continuity method). Fix any m = 1,...,M and let 

(p™,w™) e Qh x Vh satisfy the continuity method (3.1). Then (g™,u™) also 
satisfies the renormalized formulation 

f B{g^ h dx 
Jn 



-AtJ2 [ (W'^ + ^T)W'i)Wr dx 

+ At [ b{et) div fc uMh dx+ f b"(z(qZ, q™- 1 )) {el 1 - 1 } 2 <p h dx 

Jn Jn (4.1) 

rev'/? 

B"(?(g m , ($)) Ignl (<t>h)+« ■ y)l dS(x) 

= I B{g™~ 1 )(j)h dx, V0 fc G Qfc(fi), 
Jn 

for any B G C[0,oo) n C 2 (0,oo) with B(0) = and b(g) := gB'(g) - B(g). Given 
two positive real numbers a\ and a2, we use £(ai,a 2 ) and £ r (fti,a 2 ) to denote 
two corresponding numbers between a\ and a 2 that arise from second order Taylor 
expansions utilized in the proof. 

Proof. Recall the definition of m™, cf. (3.4). By taking B' '(g™)<ph as test function 
in (3.7) and repeating the proof of Lemma 5.1 in [11], we obtain (4.1) with it™ 
replaced by it m . In view of (3.5) and (3.6), this is identical to (4.1). □ 

Lemma 4.2 (Stability). Let {(gh, u h)}h>o ^ e a se Q uence of numerical solutions 
constructed according to (3.8) and Definition 3.1. For g > 0, set P(g) := ^z^g 7 . 
For any m = 1, . . . , M, we have 

/ P(e%) dx+V- YAt\\u k h \\ 2 Vh{Q) +Ar% ffusion 
Jn z 

m (4-2) 

< / P(g ) dx + Cj2^\\fh\\l H n), 
Jn fe=i 

where the numerical diffusion term M^g usion > takes the form 

^ffusion = J2 / P"(^le k h - 1 ))h k h - 1 f dx 
k=l jQ 

+ EE At [ P"(e^U] 2 T {(u k h -u)+-(u k h .^) dS(x). 

fe=irer^ r 
In particular, g h e b L°°{0, T; L 1 (0.)). 
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Proof. Since P'(p)p — P{p) = p(p) and Qh > 0, taking <\>h = 1 in (4.1) yields 

/ p( Q k h ) dx + At [ p( e k h )dw h u k h dx+ I P"{a&<t 1 ))let 1 fdx 

Jn Jn Jn 

+ At J2 f r (fl?, 9 m )) [<£] I « ■ u)+ (4.3) 

-P"(f(Q™,^))l e k h }l(u k h -^ dS(x)= [ P(g k ~ 1 )dx. 

Jn 

For k — 1, . . . , M and x € Urer 7 se * 



max{g'l(x),g k _(x)}, 1 < 7 < 2, 
mm{Q l l(x),g k _{x)}, 7 > 2, 



Q*(x) ■= 
and note that 



- P"(C r (^, #)) (Uh ■ v)l dS(x) (4.4) 

> At £ / p"( e f) [^ ((«* • „)+ - • ^) dS^). 

Next, using = u k as test function in (3.2), we obtain the estimate 

/ p(^)divfcujtdx = (^ + A)||divfc«£||i a(n) + /i||curl h u^|i, 2(n) - / f£u k h dx 
Jn Jn 

>Mll^fv,- / ft< dx, k = l,...,M. 
Jn 

(4.5) 

Applying (4.5) and (4.4) to (4.3) leads to the bound 

/ P(g k h )dx + pAt\\u k h \\ 2 Vhi n ) + f P"(Z(qI Qf)) ht't dx 
Jn Jn 

< [/(el 1 ) dx+ Yp M j n lfkh]2dx+ 2 At J n |M ^ 2 dx - 

Summing over k — 1, . . . , M yields (4.2). □ 

Since the stability estimate only provides the bound p(gh) &b L°°(Q, T; L 1 (fi)), it 
is not clear that p(gh) converges weakly to an integrable function. Hence, we shall 
next establish that the pressure is in fact uniformly bounded in L 2 (0,T; L 2 (f2)). 

To increase the readability we introduce the notation 



dx. 



Lemma 4.3 (Higher integrability on the pressure). Let {(Qh, u h)}h>o be a sequence 
of numerical solutions constructed according to (3.8) and Definition 3.1. Then 

p(Qh) e b L 2 ((0,T)xn). 
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Proof. For m = 1, . . . , M, define v™ e Vh{Ct) by 

/nJ > 



where the operator ,8 [•] is defined in (2.5). 
Since div nj^ = LI^ div, we have the identity 

div, < - ^ divB [p(^) - (p(<O n ] = - (p(Qh))a ■ 

By using t>™ as test function in the velocity method (3.2) and applying the 
previous identity, we obtain the relation 

/ P (gT) 2 dx = ( P (gT))n + (A + m) / div, < div, v%dx 
+ [ H curl, < curl, < - 

+ M E ^ / K ' "Jr K ■ "It + K x ^lr W * "lr dS ( x )- 

Repeated applications of Holder's inequality yields 

l 2 (o) 

< C (p(qT))1 + (A + M)ll div, <|| L2(n) || div, 

+ /i 1 1 curl, 1 1 £,2 (0) 1 1 curl, ug 1 1 1 L 2 (n) + 1 1 1 1 i2 (n) 1 1 < l 1 1 L 2 (n) 



E"-'/ r " 



< • f] r K* • "lr + K" x vj T K x uj T dS(x) 



rer* 

To bound the jump terms we apply Lemma 2.14: 

L 2 (f2) 



< c 



< c 



(«)>n 



+ ((1 + h5)||<|| Vh(n) + ||/ril^(o)) \\BWh) (p(eT))n] Ww^n) 
(p(Qh))1 + ((1 + hi)\\K\\v h (n) + WfTh'm) (1 + \\p{<ff)\\v>M ~ 



where the last inequality follows thanks to the estimate ||i3[</>]|| w 1 . 2 (n) < C|| < / , IU 2 (o)- 
Finally an application of Cauchy's inequality (with e) yields 

\\p(QT)\\h { a) < C (l + (p(^))o + IKIIk(o) + ll/rili 2 ( f2) ) • 

Finally we multiply this inequality by At, sum over m = 1,...,M, and apply 
Lemma 4.2. This concludes the proof. □ 

5. Convergence 

Let {(qh, u h)} h>0 be a sequence of numerical solutions constructed according 
to (3.8) and Definition 3.1. In this section we will prove that a subsequence of 
this sequence converges to a weak solution of the semi-stationary Stokes system, 
thereby proving Theorem 3.3. 

In view of Section 4, we have the following ^.-independent bounds: 

Q h G b L°°(0, T; IT>(Sl)) n L 2 \{Q, T) x fi), 
u, e b L 2 (0,T; L 2 (ft)), 
div,u, e b L 2 (0,T; L 2 (ft)), 
curl,M, G b i 2 (0,T;i 2 (r!)). 
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Consequently, we can assume that there exist limit functions g € L°°(0, T; L 7 (Q)) n 
L 2 ~t((0,T) x Q) and u G L 2 (0,T; W 1,2 (ft)) such that 

^ e in L°°(0,T;LT(n)) fl^((0,T) x SI), 

ti- 
lth - 

and, by Lemma 2.11, 



u h h ^°u in L 2 (0,T;L 2 (n)), 



(5.1) 



(5.2) 



divfttt,, ^ U divM in L 2 (0,T;L 2 (ft)), 
cwl h u h ^°curlw in L 2 (0,T;i 2 (Q)). 
Moreover, 

Q h - 1 £> 7 , 0ft £ 7+1 , Qhlogg h glogg, 

where each signifies weak convergence in a suitable L p space with p > 1 . 

Finally, ^ log p h converge respectively to g, glog g in C([0, T]; L^ eak (0)) for 
some 1 < p < 7, cf. Lemma 2.2 and also [5, 12]. In particular, g, glog g, and glog g 
belong to C([0,T];^ eak (Q)). 

5.1. Density method. 

Lemma 5.1 (Convergence of gu). Given (5.1) and (5.2), 

ghUh h ^-° gu in the sense of distributions on (0,T) x SI. 
Proof. Denote by u~h the function 

u h (t,-)=uf, tG(t m -\t m ], m=l,...,M, 

where it™ is defined in (3.4). By standard properties of the Nedelec interpolation 
operator, ||5ft||z,2( ,T;£2(Q)) < C\\u h \\ L 2 {0 ^ T . L 2 (n}} . This, (3.6), and Lemma 4.2 
allow us to conclude that u h G b L 2 (0,T; W div ' 2 (fl)) and 

m=i re r^ Jr 

Using these bounds, we can apply to (3.7) the calculations leading to Lemma 
5.6 in [11], resulting in the bound 

a t fc ( efc )ebi 1 (o,T;W- 1 - 1 (fi)). (5.3) 
At the same time, Lemma 2.13 tells us that 

\\u h (t,x) - u h (t,x - OIU 2 (o,T;£2(n { )) ^0 as |£| -> 0, uniformly in /i. (5.4) 
In view of (5.3) and (5.4), an application of Lemma 2.3 concludes the proof. □ 

Lemma 5.2 (Continuity equation). The limit pair (g,u) constructed in (5.1) and 
(5.2) is a weak solution of the continuity equation (1.1) in the sense of Definition 
24- 

Proof. Denote by the function 

u h (t,-)=uf, te(t m -\t m ], m=l,...,M, 

where ti™ is defined in (3.4). 

Fix a test function <f> e Cg°([0,T) x fi) and introduce the piecewise constant 
projections <f> h := Iljfy, 0™ := II^ m , and 0™ := 0(t, <**■ 
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By using <j>™ as test function in (3.7) and preforming the same calculations as in 
the proof of Lemma 6.4 in [11], we work out the identity 

f [ d } t l {g h )cf) h dxdt= [ [ g h UhV(t> dxdt + w{h), (5.5) 
Jo Jn Jo Jn 

where \u>(h) \ < Ch^ and 

/ / d£(Qh)<f>h dxdt 1 ^ - / / g4> t dxdt- I Q O <f>(0,x) dx, 
Jo Jn Jo Jn Jn 

where we have relied on (5.1) and the strong convergence g Q h — > g° a.e. in ft. 
Next, 

/ / Q h UhV 4> dxdt = j \ Q h u h V(j)+ Qh{uh~ u h )^4> dxdt. 
Jo Jn Jo Jn 



In view of Lemma 5.1, 

rT 







QhUh'V'^ dxdt > / / guVcj) dxdt. 



By a standard error estimate for n^ 7 (cf. [13]), 

/ / g h (ut - u h ) V</> dxdt 
Jo Jn 



< hC\\Q h \\ L 2( 0iT . L 2( n ^ || Vh w/i||L 2 (o,T;i, 2 (n)) II V0||ioo( O)T;i oo( n )) < Ch 4 , 

where the final inequality follows from Lemmas 2.12, 4.2, and 4.3. 

Summarizing, sending h — > in (5.5) delivers the desired result (2.2). □ 

5.2. Strong convergence of density approximations. To establish the strong 
convergence of the density approximations gh, we will utilize a weak continuity 
property of the effective viscous flux: P e s(Qh, Uh) = p(&h) — (A + /i) div Uh- 

To derive this property we exploit the div-curl structure of the velocity scheme 
(3.2) combined with the commutative properties (2.9) of Vh- More specifically, in 
view of the commutative property (2.9), the function Vh — 11^ VA -1 ^ satisfies 
div/j Vh — gh and curU Vh = on elements away from the boundary. The crucial 
point is that the curl part of the velocity method (3.2) vanishes when this Vh is 
utilized as a test function. 

Lemma 5.3 (Discrete effective viscous flux). Given the weak convergences listed 
in (5.1) and (5.2), 

lim / / P cS (g h ,u h ) gh<ptp dxds = / / P cS (g,u) g (f>ip dxds, 
h ^°Jo Jn Jo Jn 

for all e (ft) and V> € C°°(0,T). 

Proof. Fix (f> e Cq 30 (ft), ip £ C°°(0,T), and for each h > introduce the test 
function 

v h (-,t) = ijUX[ ( f>A[Qh-Q](;t)], te(0,T). 

where the operator A [•] is defined in (2.6). 

By virtue of (2.9) and curl .4 [•] = 0, we have the identities 

div,, v h = iPTL® (V<M [Qh - Q]) + ^h (<KQh - Q)) 

and 

curU v h = (V<f> xA[g h -g])- 
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For m = 1, . . . , M, set v™ :— ^ //m-i s ) Taking as test function 
in the velocity method (3.2), utilizing the above identities, multiplying by At, and 
summing over m, we obtain 

/ / Pes(Qh,u h )(g h - g)H> dxds 
Jo Jn 

= - / P e s(Qh, u h )V<t> ■ A [g h - g] ip + fh fiXM [g h - g]) ip dxds 
Jo Jn 

+ / yucurl/j u h (V<f> x A [qh - £»]) ip dxds 
Jo Jn 

+ n V h e ^ / ip hu h - uj r \v h ■ v\ T + \u h x i/] r K x Hlr dS(x)ds. 
rer; 7o ^ 



(5.6) 

In view of (5.3), the following h independent bounds are immediate: 
$A[Q h ] = A [d' t l g h ] G b ^(O.r;^- 1 ' 1 ^)), 
A[g h }e h L 2 (0,T;W^(fl)). 

Consequently, Lemma 2.3 can be applied with the result that (A [gh]) 2 — ^° (A [g]) 2 . 
Thus, 

A[g h -g] h ^°0, mL 2 (0,T;L 2 (n)). (5.8) 
Now, using (5.8) together with (5.1) and (5.2), we send h — > in (5.6) to obtain 

lim / / P cS (g h ,u h )(g h - g)4>ip dxds 
h ^°Jo Jn 

= ! im n ^ E he ~' I ^ I K • "] K(<M - e]) • (5 - 9) 

fc-o rer , Jo Jt 

+ K x v\ {TlX(M [g h -g])xiy] dS{x)dt. 

Lemma 2.14 yields 



fi 

rer 



V h*- 1 f <P f \u h - v\ T \U V h{4>A [g h -q])-i>} t dS(x)ds 

V h f - 1 f V / \u h x V \ T {uX(ct>A[Qh - g]) x uj r dS(x)ds 
tri Jo Jr 



fi 

< h%C\\^ L <~ {0iT) \\u h \\ L 2 {0iT . Vh{n)) \\V{(t>A[Qh - Q])\\L2(o,T;Li(n)) < Chi, 

where the last inequality follows from (5.7) and Lemmas 4.2 and 4.3. Applying the 
previous bound to (5.9) yields the desired result. □ 

We can now infer the strong convergence of the density approximations. 

Lemma 5.4 (Strong convergence of Qh)- Suppose that (5.1)-(5.2) holds. Then, 
passing to a subsequence as h — > if necessary, 

Qh — >• Q a - e - in (O,? 1 ) x ft. 

Proof. In view of Lemma 5.2, the limit (g, u) is a weak solution of the continuity 
equation and hence, by Lemma 2.6, also a renormalized solution: 

(glogg) t + div ((f?log g) u) = gdivu in the weak sense on [0,T) x ft. 
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Since t i — ► q log g is continuous with values in an appropriate Lebesgue space 
equipped with the weak topology, we can use this equation to obtain for any t > 

/ (f?l°g O) (t) dx — / g log g dx = — / / g div u dxds (5.10) 
Jn Jn Jo Jn 

Next, we specify (f>h = 1 as test function in the renormalized scheme (4.1), 

multiply by At, and sum the result over m. Making use of the convexity of z log z, 

we conclude that for any m = 1 , . . . , M 

f g% log g%dx- f g° h log g° h dx<-jTAt f g% div uf dxdt. (5.11) 
Jn Jn fc=1 Jn 

In view of the convergences stated at the beginning of this section and strong 
convergence of the initial data, we can send h — > in (5.11) to obtain 

/ ( plogf?) (t) dx — / go log g dx < — / / gdivu dxds. (5-12) 
Jn ^ ' Jn Jo Jn 

Subtracting (5.10) from (5.12) gives 

/ ( f?l°g Q ^ £"l°g Q) (t) dx < — / / g div u — g div u dxds, 
Jn^ ' Jo Jn 

for any t € (0,T). Lemma 5.3 tells us that 

J J (gdivu - gdivu) <fi dxds = ° - J J - ^pQ^j 4> dxds > 0, 

for all <p € Co°(n) n {</> > 0}, where the last inequality follows as in [5, 12], so the 
following relation holds: 



g\ogg=g\ogg a.e. in (0,T) x fi. 
Now an application of Lemma 2.1 finishes the proof. □ 
5.3. Velocity method. 

Lemma 5.5 (Velocity equation). The limit pair (g,u) constructed in (5.1)-(5.2) 
is a weak solution to the velocity equation (1.2) in the sense of Definition 2-4- 

Proof. Fix v € L 2 (0,T; w£' 2 (Sl)), and set v h = and v% = ^ £-i v h dt. 

Then, setting v™ as test function in the velocity method (3.2), multiplying with 
At, and summing over all m = 1, . . . , M, leads to the identity 

/ / (i curl h u h curl v + \) div h u h - p(g h )] divv dxdt 

Jo Jn 

rerl Jo Jr (5.13) 
+ [u fc xi/] r [(n^«) xi/] r dS(x)dt 

/ fh^Xv dxdt, 
o Jn 

where we have also used (2.9). From Lemma 5.4 and (5.1), we have that p(Qh) 
p(g) in L 2 (0, T; L 2 (£l)). Furthermore, Lemma 2.14 tells us that the jump terms 
converge to zero. Hence, we can send h — > in (5.13) to obtain that the limit (g, u) 
constructed in (5.1)-(5.2) satisfies (2.3) for all test functions v e L 2 (0, T; W 1,2 (fi)). 

□ 
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